#### relevant code from lecture 5 #####

rikz <- read.table( 'ecol 562/RIKZ.txt', header=TRUE)
rikz$richness <- apply(rikz[,2:76], 1, function(x) sum(x>0))
mod3 <- lm(richness~NAP*factor(week), data=rikz)

mod3a <- lm(richness~factor(week)+factor(week):NAP - 1, data=rikz)

out.mod <- data.frame(est=coef(mod3a), confint(mod3a), week=rep(1:4,2), 
    parm=rep(c('Intercept', 'Slope'), c(4,4)))
names(out.mod)[2:3] <- c('low95','up95')

myprepanel.ci <- function(x,y,lx,ux,subscripts,...) {
list(xlim=range(x, ux[subscripts], lx[subscripts], finite=TRUE))
}

library(lattice)
dotplot(week~est|parm, data=out.mod, xlab='Estimate', ylab='Week', 
prepanel=myprepanel.ci, lx=out.mod$low95, ux=out.mod$up95, layout=c(2,1), 
panel=function(x, y, subscripts){
panel.dotplot(x, y, col=4, pch='+', cex=.6)
panel.segments(out.mod$low95[subscripts], y, out.mod$up95[subscripts], y, 
   lwd=3, col='dodgerblue4', lineend=1)
panel.points(x, y, col='white', pch=16, cex=1.1)
panel.points(x, y, col='dodgerblue4', pch=1, cex=1.2)
panel.abline(v=0, col=2, lty=2)
}, scales=list(x='free'))

########### new code from lecture 6 ###############

#interval overlap function 
nor.func1 <- function(alpha, model, sig) 
1 - pt(-abs(qt(alpha/2,model$df.residual))*
sum(sqrt(diag(sig))) / sqrt(c(1,-1) %*% sig %*%c (1,-1)), model$df.residual) - pt(abs(qt(alpha/2,model$df.residual)) * sum(sqrt(diag(sig))) / sqrt(c(1,-1) %*% sig %*% c(1,-1)), model$df.residual, lower.tail=F)

#function whose roots we need
nor.func2 <- function(a,model,sigma) nor.func1(a,model,sigma)-.95

#obtain confidence level for intercepts
#initialize storage vector
xvec1b <- numeric(6)
#create sigma for intercepts
vcov(mod3a)[1:4,1:4]->vmat
#running index
ind <- 1
for(i in 1:3) {
  for(j in (i+1):4){
     sig <- vmat[c(i,j),c(i,j)]
     #solve for alpha
     xvec1b[ind] <- uniroot(function(x) nor.func2(x,mod3a, sig), c(.001,.999))$root
     ind<-ind+1
}}
#average the confidence levels
ci1b <- round(mean(1-xvec1b),3)

#obtain confidence level for slopes
#initialize storage
xvec1c <- numeric(6)
#create sigma for slopes
vmat <- vcov(mod3a)[5:8,5:8]
#location in storage vector
ind <- 1
for(i in 1:3) {
 for(j in (i+1):4) {
     sig <- vmat[c(i,j),c(i,j)]
     #solve for alpha
     xvec1c[ind] <- uniroot(function(x) nor.func2(x,mod3a, sig), c(.001,.999))$root
     ind <- ind+1
}}
#average the confidence levels
ci1c <- round(mean(1-xvec1c),3)

#calculate confidence intervals for slopes and intercepts
slope.ci <- confint(mod3a,level=ci1c)[5:8,]
int.ci <- confint(mod3a,level=ci1b)[1:4,]
# add CI to the data frame from before
out.mod$low50 <- c(int.ci[,1],slope.ci[,1])
out.mod$up50 <- c(int.ci[,2],slope.ci[,2])

#add a new line to the panel function from before
dotplot(week~est|parm, data=out.mod, xlab='Estimate', ylab='Week', 
prepanel=myprepanel.ci, lx=out.mod$low95, ux=out.mod$up95, layout=c(2,1), 
panel=function(x, y, subscripts){
panel.dotplot(x, y, col=4, pch='+', cex=.6)
panel.segments(out.mod$low95[subscripts], y, out.mod$up95[subscripts], y, lwd=3, col='dodgerblue4', lineend=1)
#new line
panel.segments(out.mod$low50[subscripts], y, out.mod$up50[subscripts], y, lwd=6, col='dodgerblue1', lineend=1)
panel.points(x, y, col='white', pch=16, cex=1.1)
panel.points(x, y, col='dodgerblue4', pch=1, cex=1.2)
panel.abline(v=0, col=2, lty=2)
}, scales=list(x='free'))